#
# IrtScaleVotes.R
#
# Replication file for
# Hill, Seth J. and Gregory A. Huber. 2018. ``On The Meaning of Survey Reports of Roll Call `Votes','' American Journal of Political Science.
#

rm(list=ls())
if (.Platform$OS == "windows") {
  # Set working directory in location of script.
  .doit <- function() { # only works with R.exe; trap errors if using Rscript.exe
  frame_files <- lapply(sys.frames(), function(x) x$ofile);frame_files <- Filter(Negate(is.null), frame_files)
  PATH <- dirname(frame_files[[length(frame_files)]]); setwd(PATH) ; rm(PATH,frame_files)
  }
  try(.doit(),silent=T)
}

require(foreign)
require(pscl)
require(data.table)

desired.samples <- 10e03
burnin <- 150e3
thin <- 20
maxiter <- thin*desired.samples + burnin

# Functions.
.checkConvergence <- function(res) {
  # Check convergence of an mcmc array.
  require(coda)
  cat("\nGeweke.\n")
  gd <- geweke.diag(res)
  cat("\nCumulative distribution of absolulte Geweke stats. Want ~95% less than 1.96.\n")
  print(round(cumsum(prop.table(table(round(abs(gd[[1]]),1)))),3))
  cat(sprintf("Percent less than 1.96: %1.2f.\n\n",100*mean(abs(gd[[1]]) < 1.96,na.rm=T)))

  #cat("\nRaftery-Lewis.\n")
  #print(rd <- raftery.diag(res))
}

.extractId  <- function(id) {
  # Extract posterior ideal points and make a histogram
  # from an ideal() result object.
  id.points <- data.table(id$xbar)
  setnames(id.points,c("ideal.point"))
  id.points[,V1 := rownames(id$xbar)]
  id.points[,hist(ideal.point,breaks=20,main="",xlim=c(-3,3))]
  return(id.points)
}

runIdeal <- function(dat,vote.names,normalize=T,leg.cons=NULL,...) {
  # Run ideal() irt model.
  # Arguments.
  #  dat -- data frame of votes (columns), legislators (rows)
  #  vote.names -- which column names/numbers contain votes
  #  normalize -- normalize argument passed to ideal().
  #  leg.cons -- legislator constraints when not NULL; argument
  #   `x' to pscl:::constrain.legis().
  # Returns list.
  
  # Create rollcall class object.
  rc.mat <- rollcall(dat[,vote.names],legis.names=dat$V1,vote.names=vote.names)

  # Constrain legislator ideal points if passed.
  if (!is.null(leg.cons)) {
    cl <- constrain.legis(obj=rc.mat, x=leg.cons, d = 1)
    flush.console()
    # Estimate ideal with constraints.
    id <- ideal(rc.mat,d=1,priors=cl,startvals=cl,normalize=normalize,maxiter=maxiter,thin=thin,burnin=burnin,store.item=T,...)
  } else {
    # Estimate ideal without constraints.
    id <- ideal(rc.mat,d=1,normalize=normalize,maxiter=maxiter,thin=thin,burnin=burnin,store.item=T,...)
  }
  # Grab as mcmc array.
  id.mcmc <- idealToMCMC(id)
  .checkConvergence(id.mcmc)

  # Check item params.
  item.params <- data.table(id$betabar)
  item.params[,bill.id := rownames(id$betabar)]
  setkeyv(item.params,"Discrimination D1")
  print(item.params)

  # Extract ideal point posterior means and plot.
  id.points <- .extractId(id)
  return(list("id"=id,"items"=item.params,"id.points"=id.points))
}

.extractItemParams <- function(id) {
  # Extract the posterior mean and precision for
  # the item parameters from a 1D ideal() model.
  # Returns list with two vectors, each with m
  # rows for each item in id and 2 columns for the
  # discrimination and difficulty params.
  means <- cbind(apply(id$beta[,,1],2,mean),apply(id$beta[,,2],2,mean))
  vars <- cbind(apply(id$beta[,,1],2,var),apply(id$beta[,,2],2,var))
  colnames(means) <- colnames(vars) <- c("Discrim","Diff")
  return(list("bp"=means,"bpv"=vars^-1))
}

#
# Call in data.
# Use only respondents who passed screener.
#
use.screened <- T
cat("Using only respondents who passed screener...\n")
dat <- read.csv("RawData/PreppedIRTData-Screened.csv")

vote.names <- grep("bill.",names(dat),value=T)

pdf(sprintf("figures/IRTScaledVoters%s.pdf",ifelse(use.screened,"-Screened","-Unscreened")),width=6,height=8)
par(mfrow=c(4,1),mar=c(4.1,4.1,1.6,1.1))
#
# =============
# First scale all respondents together.
# =============
#
cat("Running ideal() on all survey respondents...\n");flush.console()
all.id <- runIdeal(dat[dat$assignedCondition != "Senate",],vote.names)
title(main="All respondents")

#
# =============
# Then separately by assigned condition/Senate.
# =============
#
cat("Running ideal() on standard condition only...\n");flush.console()
std.id <- runIdeal(dat[dat$assignedCondition == "Standard",],vote.names)
title(main="Respondents in standard condition")

cat("Running ideal() on informational condition only...\n");flush.console()
inf.id <- runIdeal(dat[dat$assignedCondition == "Info",],vote.names)
title(main="Respondents in informational condition")

cat("Running ideal() on Senators only...\n");flush.console()
sen.id <- runIdeal(dat[dat$assignedCondition == "Senator",],vote.names)
title(main="Senators on same 12 roll calls")

dev.off()

#
# =============
# Joint scalings.
# =============
#

cat("Projecting standard respondents into Senators-only space...\n");flush.console()

# Estimate on respondents joint with Senators; 
# fix item parameters to posterior means from
# Senate-only scaling.
# Extract posterior mean item params.
items <- .extractItemParams(sen.id$id)
# Create list of priors. Huge precision on each item
# parameter leads to ~fixed item params.
priors <- list("bp"=items$bp,"bpv"=100e3) 

# Joint scale with standard respondents fixing item
# parameters.
sen.std.id <- runIdeal(dat[dat$assignedCondition %in% c("Senator","Standard"),],vote.names,normalize=F,priors=priors)
title(main="Senators and standard condition")
# Check: should be same after estimation!
print(items$bp)

cat("Projecting info respondents into Senators-only space...\n");flush.console()
# Join scale with info respondents fixing legislator
# ideal points.
sen.inf.id <- runIdeal(dat[dat$assignedCondition %in% c("Senator","Info"),],vote.names,normalize=F,priors=priors)
title(main="Senators and informational condition")
# Check: should be same after estimation!
print(items$bp)
  
# Save results.
save(dat, vote.names, all.id, std.id, inf.id, sen.id, sen.std.id, sen.inf.id, file=sprintf("RawData/IRTModelResults%s.RData",ifelse(use.screened,"-Screened","-Unscreened")))

#
# Call in data.
# Use all respondents.
#
use.screened <- F
dat <- read.csv("RawData/PreppedIRTData-Unscreened.csv")

vote.names <- grep("bill.",names(dat),value=T)

pdf(sprintf("figures/IRTScaledVoters%s.pdf",ifelse(use.screened,"-Screened","-Unscreened")),width=6,height=8)
par(mfrow=c(4,1),mar=c(4.1,4.1,1.6,1.1))
#
# =============
# First scale all respondents together.
# =============
#
cat("Running ideal() on all survey respondents...\n");flush.console()
all.id <- runIdeal(dat[dat$assignedCondition != "Senate",],vote.names)
title(main="All respondents")

#
# =============
# Then separately by assigned condition/Senate.
# =============
#
cat("Running ideal() on standard condition only...\n");flush.console()
std.id <- runIdeal(dat[dat$assignedCondition == "Standard",],vote.names)
title(main="Respondents in standard condition")

cat("Running ideal() on informational condition only...\n");flush.console()
inf.id <- runIdeal(dat[dat$assignedCondition == "Info",],vote.names)
title(main="Respondents in informational condition")

cat("Running ideal() on Senators only...\n");flush.console()
sen.id <- runIdeal(dat[dat$assignedCondition == "Senator",],vote.names)
title(main="Senators on same 12 roll calls")

dev.off()

#
# =============
# Joint scalings.
# =============
#

cat("Projecting standard respondents into Senators-only space...\n");flush.console()

# Estimate on respondents joint with Senators; 
# fix item parameters to posterior means from
# Senate-only scaling.
# Extract posterior mean item params.
items <- .extractItemParams(sen.id$id)
# Create list of priors. Huge precision on each item
# parameter leads to ~fixed item params.
priors <- list("bp"=items$bp,"bpv"=100e3) 

# Joint scale with standard respondents fixing item
# parameters.
sen.std.id <- runIdeal(dat[dat$assignedCondition %in% c("Senator","Standard"),],vote.names,normalize=F,priors=priors)
title(main="Senators and standard condition")
# Check: should be same after estimation!
print(items$bp)

cat("Projecting info respondents into Senators-only space...\n");flush.console()
# Join scale with info respondents fixing legislator
# ideal points.
sen.inf.id <- runIdeal(dat[dat$assignedCondition %in% c("Senator","Info"),],vote.names,normalize=F,priors=priors)
title(main="Senators and informational condition")
# Check: should be same after estimation!
print(items$bp)
  
# Save results.
save(dat, vote.names, all.id, std.id, inf.id, sen.id, sen.std.id, sen.inf.id, file=sprintf("RawData/IRTModelResults%s.RData",ifelse(use.screened,"-Screened","-Unscreened")))
